Creating Artificial Ice States Using Vortices in Nanostructured Superconductors 
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We demonstrate that it is possible to realize vortex ice states that are analogous to square and 
kagome ice. With numerical simulations, we show that the system can be brought into a state that 
obeys either global or local ice rules by applying an external current according to an annealing 
protocol. We explore the breakdown of the ice rules due to disorder in the nanostructure array and 
show that in square ice, topological defects appear along grain boundaries, while in kagome ice, 
individual defects appear. We argue that the vortex system offers significant advantages over other 
artificial ice systems. 
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Geometric frustration occurs when a system is con- 
strained by geometry in such a way that the pairwise 
interaction energy cannot be simultaneously minimized 
for all constituents, and appears in water ice 1\, spin 
systems [3, [3, , an d a variety of other systems in both 
physics [3] and biology [f|. A specific example of frus- 
tration occurs in the classical spin ice system where the 
constituents of the system are magnetic spins on a grid 
of corner-sharing tetrahedra. The spins are constrained 
to point along the lines connecting the middle points of 
the tetrahedra [3, 4] and pairs of spins can minimize their 
energy by adopting a head-to-tail configuration. It is not, 
however, possible for the four spins on a tetrahedron to 
simultaneously satisfy each of the six pairwise interac- 
tions in a head-to-tail fashion; the best the system can 
do is to satisfy four interactions out of six, leaving two 
pairs in a head-to-head or tail-to-tail configuration. As 
a result, in the ground state configuration each tetrahe- 
dron obeys the so-called "ice rule" of a two-in two-out 
configuration with two spins pointing toward the center 
of the tetrahedron and two spins pointing away from it. 
Defects appear in the form of magnetic monopoles [7|. 

Recently, there has been growing interest in crea ting 
model systems that exhibit spin ice behavior [3, [3, E3, 
E1E3, M , E3] and that allow the individual constituents 
to be imaged directly, unlike molecular or atomic ices. 
For example, Wang et al. [3] created artificial square 
ice using single-domain rectangular ferromagnetic islands 
arranged in a square lattice such that four islands meet at 
every vertex point. They found that as the inter-island 
interaction increased, the system preferentially formed 
ice-rule-obeying vertices, but it did not reproduce the 
known ground state of two-dimensional (2D) spin ice, 
where the two "in" magnetic moments are on opposite 
sides of the vertex. This could be due to the relative 
weakness of the magnetic interactions. It has recently 
been shown that certain dynamical annealing protocols 
permit the system to approach the ground state more 
closely [3, E3] • Similar studies have been performed for 
a 2D kagome ice system [13, EH where the local ice-rules 
were obeyed and defects such as three-in or three-out 
were absent [13]. In the colloidal artificial ice system of 
Ref. [13], the local dynamics can be accessed easily via 



video microscopy; however, the ice arrays in this system 
are limited to relatively small sizes in experiment. 

Here we propose that a particularly promising artifi- 
cial ice system could be created using vortices in super- 
conductors with appropriately designed nanostructured 
arrays of artificial pinning sites. There has been ex- 
tensive experimental work showing that a rich variety 
of different pinning array geometries can be fabricated 
[13, E3, E3, E3, E3, H3|, and various types of experi- 
mental techniques exist for directly imaging vortices in 
these arrays [13, E3, E3, IS! • The vortex system has sev- 
eral advantages over other artificial ice systems. The 
vortex-vortex interaction strength is large, permitting 
the ground state to be reached much more readily than 
in the nanomagnetic systems. An applied external cur- 
rent permits the straightforward realization of different 
dynamical annealing protocols. New types of defects can 
be studied by merely increasing or decreasing the mag- 
netic field to create vacancies or interstitials that locally 
break the ice rules, while transport properties and criti- 
cal currents can be measured which are not accessible in 
the other systems. 

To form square vortex ice, we propose using an ar- 
rangement of elongated double- well pinning sites. Non- 
superconducting islands with the double-hump shape il- 
lustrated in Fig. 1(a) placed within a superconducting 
layer have a pair of potential minima at the highest points 
of the island where the superconducting layer is the shal- 
lowest. A single vortex trapped over each island will sit 
at one of the two minima, depending on the interactions 
with nearby vortices. By changing the arrangement of 
the islands, different types of ice can be created. For 
square ice, shown in Fig. 1(a), four islands come together 
at each vertex and the state of each island is defined as 
"in" if the vortex sits close to the vertex and "out" oth- 
erwise. We define TI[yi clS the number of "in" vortices at a 
vertex. In Fig. 1(a), the vortices have formed an n- m = 2 
ice-rule-obeying ground state configuration. Figure 1(b) 
shows a kagome spin ice arrangement with three islands 
surrounding each vertex. In this case, the lowest energy 
state has n; n = 1 or n- m — 2 at each vertex, but there is 
no overall ordering into a unique ground state. 

To study the vortex ice, we perform numerical Simula- 
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FIG. 1: Schematic of the nanostructured pinning site con- 
figurations producing ice states. Double-lobed objects: pins; 
open mesh objects: vortices, a) Square ice ground state, b) 
One possible biased ground state of the kagome ice system. 



tions of a 2D sample with periodic boundaries contain- 
ing N p elongated pinning sites in the square or kagome 
configurations illustrated in Fig. 1 and N v = N p vor- 
tices. A vortex i at position R^ obeys the following 
overdamped equation of motion: r){dHi/dt) — f? v + 
ff + i d + if . The damping constant 77 = 4>f ) d/2n( t 2 p N , 
where 0o = h/2e is the flux quantum, £ is the super- 
conducting coherence length, is the normal state re- 
sistivity of the material, and d is the thickness of the 
superconducting crystal. The vortex-vortex interaction 

force is given by f™ = Yfj^i foKi(Rij I 'ty&ij , where 
K\ is the modified Bcsscl function appropriate for stiff 
three-dimensional vortex lines, A is the London penetra- 
tion depth, /o = 4>Q/(2irpo\ 3 ), Rij — |Rj — R,j\, and 
Ry = (R; — Rj)/Rij. The substrate force f/ arises 
from the elongated pins, (? k = Y, k " fo(fp/r p )Rf k &(r p - 

Rt k )&tk + fo(U/r P )Rie(r P ~ Rk)&ik + fo(fb/W ~ 
RlMl ~ 4)RL Here R± = |R, - R* ± Jpf |, R^ 1 = 
|(Ri — Rfc) • p^ || I, R^ is the position of pin k, and p| 

(Px) is a unit vector parallel (perpendicular) to the axis 
of pin k. Each vortex is constrained to stay within a pin 
composed of two half-parabolic wells of radius r p = 0.4A 
separated by an elongated region of length 21 which con- 
fines the vortex perpendicular to the pin axis and has a 
repulsive potential or barrier of strength /& parallel to 
the axis which pushes the vortex out of the middle of the 
pin into one of the ends. We take I = 2/3A or 5/6A and 
vary the lattice constant a of the pinning array between 
a = 2. OA and 8. OA. The driving force id represents the 
Lorentz force from an applied current. The thermal force 
f 2 T comes from thermal Langevin kicks and is set to zero 
except during the annealing of the kagome ice. 

Square Ice - We prepare the square ice system using 
a dynamical annealing procedure inspired by the nano- 
magnetic ice results of Refs. 0, In our simulations, 
we place one vortex in each pin at a random position and 
then use a protocol of a rotating in-plane applied current 
with decreasing amplitude, i d = A ac (t)(cos(2Trt/T r )x + 
sin(27rt/T r )y), where T r — 1000 simulation time steps, 
A ac (t) = ±(A - SA[t/5t\), A = 2.0/o, St = 10000 sim- 
ulation time steps, and SA = 0.0f/o. The force direction 
is reversed each time the magnitude of the force is de- 



creased. We measure the number of vertices of each type 
that appear after completing the dynamical annealing. 
For the kagome ice system, we obtain the vortex config- 
urations from standard thermal simulated annealing. 

To determine how effectively the dynamical anneal- 
ing protocol brings the square ice system to the ground 
state, we introduce disorder to the system by replacing 
the delta-function distributed barriers fb at the center of 
each pinning site with barriers of normally distributed 
strength, where the mean strength is /& and the width 
of the distribution is a. In Fig. 2 we illustrate the ver- 
tices that have reached the ground state configuration of 
7ii n = 2 in a square ice sample with a = 2.5A, I = 5/6A 
and fb = 0.25/o for differing disorder widths a. The dots 
represent vertices in the ice-rule obeying ground state, 
while the closed black circles indicate higher energy ver- 
tices that we term ice-rule defects Dj since they still obey 
the n; n = 2 ice rule but have the two "in" vortices ad- 
jacent to one another. The open circles mark the high- 
est energy vertices that we term non-ice-rule defects Dni 
since they do not obey the n- m — 2 ice rule but have, for 
example, n- m = 3 or nj n = 0. For a < 0.1, the system can 
reach the ordered ground state as shown in Fig. 2(a). As 
the central barriers of the pins become more nonuniform 
with increasing a, some pinning centers act as nucleation 
sites for grain boundaries, as illustrated in Figs. 2(b,c) 
for (7 = 0.1 and a = 0.5. In general, we find that for 
O.f < a < 0.7, all of the defected vertices form closed loop 
grain boundaries and the ratio of D\ to Dni is f:l due 
to geometric constraints. For a > 0.7, Fig. 2(d) shows 
that a proliferation of £>ni occurs so that the Dni out- 
number the D\. The grain boundary loops interact and 
wind around the sample, making it difficult to determine 
the relation between a and the grain boundary length. 
We find that individual Z?ni can appear outside of grain 
boundaries, while D\ always remain confined to grain 
boundaries, suggesting that there could be a disorder- 
induced phase transition when the Dni proliferate. We 
also find that doubly occupied pinning sites with two vor- 
tices each can act as grain boundary nucleation sites, as 
illustrated in the inset of Fig. 4(b). 

In Figure 3(a), we plot the percentage of vertices Pes 
that are in the ice-rule obeying ground state as a func- 
tion of time during the dynamical annealing procedure 
in a sample with a = 2.5A, I — 5/6A, fb = 0.25/o, 
and different values of a. At early times, when |A oc | 
is close to j4o> a ll of the vortices follow the drive and 
switch back and forth inside the pinning sites. As \A ac \ 
decreases, a transition occurs when the vortices cease to 
follow the driving direction and become locked into one 
position in the pinning site. For a = 0, this locking tran- 
sition is relatively sharp and occurs at \A ac \ as 0.82/o. 
Nonzero values of a broaden the transition significantly 
and cause some vertices to lock into the ground state 
at much earlier times; at the same time, complete lock- 
ing of all vertices into the ground state can no longer be 
achieved within the finite time of the dynamical anneal- 
ing process. We quantify the broadening of the transition 
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FIG. 2: Grain boundary images in square ice samples with 
a = 2.5A, I = 5/6A, and ft = 0.25 for increasing disorder 
width a. Dots: ground state m n = 2 ice-rule obeying vertices; 
filled black circles: ice-rule defects D\\ white circles: non-ice- 
rule defects D m . (a) a = 0. (b) a = 0.1. (c) a = 0.5. (d) 
a = 1.0. 



with increasing a by fitting the curves in Fig. 3(a) to the 
form Pcs{t) = 1 — exp(£/r). Figure 3(b) shows the fit- 
ted relaxation time r as a function of a and indicates 
the occurrence of an increasingly slow locking process as 
the disorder width increases. The dependence of Pes on 
both a and a is summarized in Fig. 3(d) for a system 
with fb — 1.0 and I = 2/3 A. Here, Pqs decreases both 
with increasing a and with increasing a as the relative 
strength of the vortex- vortex interactions decreases. 

Depending on the system parameters, it is not al- 
ways necessary to perform a dynamical annealing pro- 
cedure in order to reach the ground state. To demon- 
strate this, we prepare the sample in a random state 
and then apply a fixed amplitude rotating drive, i d — 
A(cos(2TTt/T r )x + sin(27rf/T r )y), with A = 0.01/ and 
T r = 1000 simulation time steps, for 2 x 10 6 simulation 
time steps. When the central barrier in the pin /& is 
weak, the system can reach the ordered ground state un- 
der the weak external shaking. For larger fb, the system 
cannot reach the ordered ground state without dynamical 
annealing. This is shown in Fig. 3(c), where we plot the 
final Pqs at the end of the simulation time versus fb for 
samples with a = 0.01 and varied pinning lattice constant 
a = 2. OA, 2.5A, and 3. OA. For large fb, the sample is im- 
mediately frozen into the disordered initial configuration, 
and Pqs ~ 0.125, consistent with the value expected in 
a completely random sample. As fb is lowered, a spon- 
taneous rearrangement into a partially ordered state be- 
comes possible and Pes > 0.125. The value of fb at 
which the spontaneous ordering appears increases with 
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FIG. 3: (a) Percentage Pes of ice-rule obeying ground state 
vertices vs time during the dynamical annealing process for 
different disorder widths a. From upper right to lower right, 
a = 0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, and 1.0. 
Here, a — 2.5A, / = 5/6A, and fb = 0.25. (b) Relaxation time 
r vs a for the same system, (c) Final value of Pas vs ft in 
samples subjected to a small shaking field with no dynamical 
annealing. Here I = 5/6A, a = 0.1, and a — 2. OA (open 
circles), 2.5A (filled squares), and 3. OA (open diamonds), (d) 
Pgs vs a and a in a sample with /;, = 1.0 and / = 2/3A. 



decreasing a, indicating that as the vortex-vortex inter- 
actions grow stronger in the denser pinning arrays, the 
ordered ground state is much more energetically favored. 

Kagome ice - The kagome lattice illustrated in 
Fig. 1(b) has a distinct set of ice rules from the square lat- 
tice. High energy vertices with n ln = or 3 are avoided in 
favor of the kagome- ice-rule obeying vertices with n- m = 1 
or 2. This system can form a non-unique ordered ground 
state, but only in the presence of an external biasing field. 
In Fig. 4(a) we show one possible biased ordered ground 
state for a kagome lattice with fb = 1.0 and a — ob- 
tained by applying a constant drive f d = 0.01/o(x + y) 
along a lattice symmetry direction while performing sim- 
ulated annealing. In the absence of the biasing force, 
some high energy defect vertices which take the form of 
monopoles appear in the system and there is no overall 
order, as illustrated in Fig. 4(c). We find that the kagome 
ice is more robust against the effects of disorder than the 
square ice, in agreement with experimental findings for 
nanomagnetic kagome ice The defect patterns are 

distinct from the square ice since no grain boundary state 
forms for the kagome ice due to the lack of an ordered 
ground state. Unlike the bipartite square lattice, the non- 
bipartite kagome lattice is not topologically constrained, 
making our system more closely resemble the ice state 
studied in Ref . [12] than that considered in Ref. [23[ . Al- 
though Fig. 4(c) shows that there is some tendency for 
the defected vertices to form pairs, there are no extended 
defect patterns of the type seen in Fig. 2. Since the ice 
rules in this system are enforced by the vortex-vortex 
interaction energies, we can weaken the enforcement of 
the ice rules by increasing the spacing a between pinning 



4 



(a) 
80 

70 

60 

50 

40 

(c) 



"■"Ho 

o@ 
_© o 

om 
© o 

o© 
-® o 

o© 

£3 O 
OS 
-© O 
O © 

© o 

OS 

-© o 
o© 
© o 

Of 

ti lo- 



1 t>©' a 
O® OQ 

© o© o 
o© o© 

© o© o 
o© o© 

® CO o 
OS OS 

o© o 
o© oa 

© o© o 
o© o© 

S OS o 
OS o© 

© o© o 
o© o© 

© o© o 

O ©I OS, 

o i d 



^51 1 o L 

CO o© 
© o © oja 
oo o o 

© OO OS 
OO OO 

© O© OS 
OS OS 

© oe> o® 
o © o® 

® om o-@ 
o© o© 

© o © o © 

OS OS 

© o© o-e 
o© o© 
© o© o© 
o© i o© | 
O l ' 0"H 



50 60 70 80 90 



120 
100 
80 
60 
40 
20 





IpO O OO, 



+ o o 

OQ _ O C 

O OO O OO o 



zvwr-Tj-w-qo-w 



. rfo ol 

OT3Q O O f 
„BO, O „ 00,0 

oo if ^o 3> 
o s> so o o 
a o o oo o 

oo 

Z a o oo o c 

5 O „ O „ OO OO, C 



o°oo°V* c 

i _ OO OO 



? Oo'°p „ OC 
OO OO 



, 00„ o oo. 



o o oo o o 
OO O O _ OO, 

o bo oo o -a 

OO O O O 
P O OO OO „ ( 



3 ° J° 

o „"po ■ 



o o oo o 
c? a a 
O <* 

oo oo o 



?JD G "SO" ~Oi 



o oo^ °"Z o~ oq, 

O O Ou m 
q, o o „ 00,0 O o 

■DO OO OO 
3 O O O °3 
P° ^ > « ,P 

3 „ o O OO o o 
o„ q o o , o p 



J3 + O OOO 

* °°oo J 00 o^ oo° o o°o°+ 

O OO o o oo 
- „o oq, rt q,„ o q, n oo a 

o o 00 oo o o 
OWL o o o o oo o 
„ u q,> rt a o o oq, oq, o— 
o OO o O 

„ o „oo o_oo o o o„ o o „ 

1 oo _o *o o o oa 
00 o o o o 

o oo oo o oo a „ oq, o 
oq, oq, o p _ o o oq, o* 



„2>?„ <£ C 



3a o rt q, o oq,^ c 
, oC qo oo 

3 O OO OO OC 

MS 







50 



100 

x 



150 



200 



FIG. 4: (a) Ordered biased ground state in a sample with 
kagome pinning, /(, = 1.0, / = 2/3A, and a = 3A. Open 
circles: ni n = 1 vertices; shaded circles: ni n = 2 vertices, (b) 
Percentage Pv of each vertex type vs a. Crosses: n- m = 0; 
open circles: n ln = 1; shaded circles: n ln = 2; filled squares: 
n ln = 3. Inset: grain boundary image in square ice sample 
with two doubly occupied pins (open squares) with the same 
symbols as in Fig. 2. (c) Vertex configuration after thermal 
annealing in a sample with a = 3.5A, I — 2/3, and /(, = 1.0. 
Symbols are the same as in the main panel of (b). 



sites. Figure 4(b) shows that as a increases, the system 
passes from a limit in which only kagome-ice-rule obey- 
ing vortices appear for a < 4A to a limit a > 8A where 
the vertices assume a completely random arrangement. 



In the random limit we expect to find each of the two de- 
fect vertex types with probability 1/8 and each of the two 
kagome-ice-rule obeying vertices with probability 3/8. 

There are other arrays that would obey ice-rule type 
constraints; however, the simplest cases for 2D are the 
square and kagome arrays. Previous studies of super- 
conducting wire networks arranged in kagome configura- 
tions found geometrical frustration which produced dis- 
ordered ground states [HI; however, such a system does 
not specifically have ice-rule obeying states. The artificial 
ice vortex system proposed here can be used to study the 
effect of ice-rule and non-ice- rule configurations on trans- 
port and magnetization properties, and it would also be 
possible to examine higher matching fields to see whether 
new types of ordered or disordered states appear. 

In summary, we propose that square and kagome vor- 
tex ice can be realized in nanostructured superconduc- 
tors. By using an annealing protocol of a rotating exter- 
nally applied current, the system can reach or approach 
the square ice ground state. In the presence of quenched 
disorder, defects appear in an ordered ground state back- 
ground. For moderate disorder in the square ice system, 
all of the defects are bound to grain boundaries, while for 
strong disorder, individual high energy vertices prolifer- 
ate. For kagome ice, we find no grain boundary phase 
in the presence of disorder. We predict that if the bar- 
rier for vortex motion across the center of each artificial 
pinning site is weak, the system will spontaneously orga- 
nize into a partially ordered state even without use of an 
annealing protocol. This system could have interesting 
transport and memory effects which may manifest them- 
selves as changes in the critical current, an effect which 
cannot be accessed readily in other artificial ice systems. 

We thank C. Nisoli for a useful discussion. This work 
was carried out under the NNSA of the U.S. DoE at 
LANL under Contract No. DE-AC52-06NA25396. 



[1] L. Pauling, J. Am. Chem. Soc. 57, 2680 (1935). 
[2] P.W. Anderson, Phys. Rev. 102, 1008 (1956). 
[3] R. Moessner and A. P. Ramirez, Phys. Today, 59(2), 24 
(2006). 

[4] A.P. Ramirez et al, Nature (Lo ndon) 399, 333 (1999). 
[5] Y. Han et al, larXiv:0807.3905l 

[6] H. Frauenfelder, P.G. Wolynes, and R.H. Austin, 
Rev. Mod. Phys. 71, S419 (1999); H.M. Harreis, 
C.N. Likos, and H. Lowen, Biophys. J. 84, 3607 (2003). 
[7] C. Castelnovo, R. Moessner, and S.L. Sondhi, Nature 

451. 42 (2008): L.A.S. Mol et al.. larXiv:0809.2105l 
[8] R.F. Wang et al, Nature (London) 439, 303 (2006). 
[9] C. Nisoli et al, Phys. Rev. Lett. 98, 217203 (2007). 
[10] X. Ke et al, Phys. Rev. Lett. 101, 037205 (2008). 
[11] G. Moller and R. Moessner, Phys. Rev. Lett. 96, 237202 
(2006). 

[12] M. Tanaka et al, Phys. Rev. B 73, 052411 (2006). 

[13] Y. Qi, T. Brintlinger, and J. Cumings, Phys. Rev. B 77, 



094418 (2008). 

[14] A. Libal, C. Reichhardt, and C.J. Olson Reichhardt, 

Phys. Rev. Lett. 97, 228302 (2006). 
[15] M. Baert et al, Phys. Rev. Lett. 74, 3269 (1995). 
[16] J.I. Martm et al, Phys. Rev. Lett. 83, 1022 (1999). 
[17] K. Harada et al, Science 274, 1167 (1996). 
[18] S.B. Field et al, Phys. Rev. Lett. 88, 067003 (2002); 

A.N. Grigorenko et al, ibid. 90, 237001 (2003). 
[19] G. Karapetrov et al, Phys. Rev. Lett. 95, 167002 (2005). 
[20] G. Karapetrov et al, Appl. Phys. Lett. 87, 162515 (2005). 
[21] I.V. Grigorieva et al, Phys. Rev. Lett. 99, 147003 (2007). 
[22] A.S. Wills, R. Ballou, and C. Lacroix, Phys. Rev. B 66, 

144407 (2002). 

[23] R. Moessner and S.L. Sondhi, Phys. Rev. B 68, 064411 
(2003); S.V. Isakov, R. Moessner, and S.L. Sondhi, 
Phys. Rev. Lett. 95, 217201 (2005). 

[24] M.S. Rzchowski, Phys. Rev. B 55, 11745 (1997); D. Davi- 
dovic et al, ibid. 55, 6518 (1997); K. Park and D.A. 



■5 



Huse, ibid. 64, 134522 (2001); M.J. Higgins et al, ibid. 
61, R894 (2000); Y. Xiao et al, ibid. 65, 214503 (2002). 



